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ABSTRACT 


Seismie  monitoring  requires  aeeurate  souree  characterization  in  real  time.  Aeeurate  3D  earth  models  are  essential  for 
accurate  predictions  of  seismie  observ  ables  and  source  characterization.  While  recovering  the  true  earth  structure 
has  always  been  the  goal  of  tomographic  inversions,  methods  based  on  full  wavefields  in  3D  reference  models  have 
several  distinet  advantages.  Made  possible  by  advances  in  high-performance  computation,  this  full-wave  approach 
aeeounts  for  complex  wave  propagation  in  3D  heterogeneous  earth,  enables  fuller  utilization  of  seismie  records,  and 
allows  us  to  linearize  the  inverse  problem  by  iteratively  updating  the  3D  reference  model.  An  important  benefit  of 
physieally  realistic  and  accurate  modeling  of  wavefields  in  3D  models  is  the  consistency  of  the  system  of  equations 
in  inversion.  This  is  particularly  important  for  the  integration  of  different  types  of  observations  (P,  S,  and  surfaee 
waves  including  empirical  Green’s  functions  derived  from  ambient  noise)  in  inversion. 

Two  methods  have  been  developed  in  this  field  to  carry  out  full-wave  tomography  iteratively  with  3D  reference 
models.  One  is  the  adjoint-wavefield  (AW)  method,  which  baek-propagates  the  data  from  the  receivers  to  image  the 
structure.  The  other  is  the  seattering-integral  (SI)  method,  whieh  constructs  the  strain  Green  tensor  (SGT)  databases 
and  calculates  sensitivity  kernels  for  eaeh  data  functional.  Both  methods  are  based  on  the  full-wave  theory,  and  the 
main  differences  are  in  the  computational  approaches.  In  general,  the  adjoint  method  is  CPU  intensive,  while  the 
SGT-database  approach  requires  a  large  disk  space  and  fast  network  To  date,  the  AW  and  SI  methods  use  two 
different  approaches  in  calculating  synthetic  waveforms,  i.e.  the  speetral-element  method  and  finite-difference  (FD) 
method,  respectively.  One  notable  advantage  of  the  SI  method  is  that  the  SGT  databases  make  it  possible  to  use  3D 
synthetics  in  real  time  seismie  monitoring. 

We  have  systematically  developed  a  unified  source  mechanism  and  tomographic  inversion  method  based  on  3D 
finite-difference  SGT  (FDSGT)  databases.  The  main  developments  include: 

•  a  finite-difference  code  in  the  spherieal  coordinate  system  for  regional  and  global  wave  propagation; 

•  a  new  and  efficient  implementation  of  the  Perfectly  Matched  Layer  that  prov  ides  a  better  performance  for 
waves  incident  at  a  near-grazing  angle,  low  frequency  waves  or  evanescent  waves; 

•  aeeurate  representation  of  Pn/Sn  finite-frequency  sensitivity  kernels; 

•  a  fuller  utilization  of  an  arrival  on  all  three  components  of  seismie  records; 

•  the  dependence  of  P  waveforms  on  5-wave  speed  anomalies  and  Rayleigh  waveforms  on  P-wave  speed 
anomalies; 

•  moment  tensor  inversion  using  SGT  databases;  and 

•  joint  P  and  S  velocity  tomography  using  an  integrated  body  and  surfaee  wave  dataset. 

We  have  applied  this  method  to  Iran  and  the  southeast  Tibetan  plateau  and  demonstrated  our  SGT-based,  unified 
tomography  and  moment  tensor  inversion.  The  FDSGT  moment  tensor  solutions  agree  well  with  the  global  eentroid 
moment  tensors  (GCMT).  Beeause  of  loeal/regional  stations,  3D  models,  and  shorter-period  waves  than  in  GCMT, 
FDSGT  moment  tensor  solutions  provide  better  locations,  particularly  in  constraining  depth.  The  joint  P  and  S 
velocity  model  of  the  southeast  Tibetan  plateau  from  tomography  of  - 10,000  finite- frequency  body  and  surface 
wave  measurements  shows  an  unprecedented  resolution  compared  to  previous  models  and  tomography  of  similar 
seales. 


20090914186 


201 


2009  Monitoring  Research  Review:  Ground-Based  Nuclear  Explosion  Monitoring  Technologies 


objectives 

The  objectives  of  this  work  are  to  obtain  Finite-Frequency  Seismic  Tomography  (FFST)  velocity  models  for  Eurasia 
and  to  refine  crustal  and  shallow  upper  mantle  velocity  and  attenuation  models  for  focused  areas  of  interest. 

RESEARCH  ACCOMPLISHED 

Accurate  3D  earth  models  arc  essential  for  accurate  predictions  of  seismic  observables  and  source  characterization. 
While  recovering  the  true  form  and  magnitude  of  earth  structures  has  always  been  the  goal  of  tomographic 
inversions,  recent  studies  show  that  methods  based  on  full  wavefields  in  3D  reference  models  have  several  distinct 
advantages  in  relating  seismic  observations  to  highly  heterogeneous  Earth  structures.  Made  possible  by  advances  in 
high-performance  computation,  this  full-wavefield  approach  accounts  for  complex  wave  propagation  in  3D 
heterogeneous  Earth  (Zhao  et  al.,  2005;  Tromp  et  al.,  2005;  Liu  and  Tromp,  2006;  Zhang  et  al.,  2007),  enables  fuller 
utilization  of  an  arrival  on  all  three  components  of  seismic  records  (Shen  et  al.,  2008a),  and  allows  us  to  linearize  the 
inverse  problem  by  iteratively  updating  the  3D  reference  model.  An  important  benefit  of  physically  realistic  and 
accurate  modeling  of  full  wavefields  in  3D  models  is  the  consistency  of  the  system  of  equations  in  inversion.  This  is 
particularly  important  for  the  integration  of  different  types  of  observations  (/\  S',  surface  waves  including  empirical 
Green’s  functions  derived  from  ambient  noise)  in  inversion,  which  is  an  essential  step  to  obtain  a  coherent  and  self- 
consistcnt  model  of  the  crust  and  mantle.  Seismic  tomography  based  on  path-  or  structure-averaging  (rays)  and  ID 
or  2D  reference  models  does  not  accurately  place  the  3D  sensitivities  of  seismic  observ  ations  in  a  complex  3D 
model.  The  mismatch  in  sensitivities  due  to  approximations  in  theory  becomes  problematic  when  the  target 
resolutions  are  at  scales  smaller  than  the  Fresnel  zone.  For  scale,  a  1-Hz  Pn  wave  recorded  by  a  receiver  -1000  km 
from  the  source  has  a  Fresnel  zone  width  of  -200  km  in  the  middle  of  the  wave  path  (Zhang  et  al.,  2007).  van  dcr 
Lee  and  Frederiksen  (2005)  recognize  that  a  more  consistent  and  linear  system  of  equations  offered  by  the  full- 
wavefield  approach  would  fit  the  combined  datasets  better  under  comparable  regularization  biases. 

Two  methods  have  been  developed  to  earry  out  full-wave  tomography  iteratively  with  3D  reference  models.  One  is 
the  AWr  method,  which  back-propagates  the  data  from  the  receivers  to  image  structure  (Tromp  et  al.,  2005;  Liu  and 
Tromp,  2006).  The  other  is  the  SGT-based,  SI  method,  which  calculates  and  stores  the  sensitivity  kernels  for  each 
data  functional  (Zhao  et  al.,  2005;  Zhang  et  al.,  2007).  Both  methods  are  based  on  the  full  wave  theory.  The  main 
differences  are  in  the  computational  approaches.  In  general,  the  adjoint  method  is  CPU  intensive,  while  the  SGT 
database  approach  requires  a  large  disk  spaee  and  a  fast  network.  A  notable  advantage  of  the  SI  method  is  that  it 
provides  both  the  Hessian  and  gradient  of  the  misfit  function  in  the  tomographic  inversion  and  requires  a  fewer 
number  of  simulations  per  model  iteration  to  achieve  the  same  level  of  data  variance  reduction  (Chen  et  al.,  2007). 
The  SGT  databases  can  also  be  used  in  source  moment  tensor  inversions,  making  it  possible  to  use  3D  synthetics  in 
real  time  seismic  monitoring.  To  date,  the  AW  and  SI  methods  use  two  different  approaches  to  calculate  synthetic 
waveforms:  the  spcctral-clcmcnt  method  (SEM)  (Komatitseh  et  al.,  2002)  and  finite-diffcrcncc  method  (Olsen, 

1904;  Zhang  ct  al.,  2008;  Zhang  and  Shen,  2009;  Zhang  ct  al.,  2009),  respectively. 

lTnified  tomography  and  source  inversions  based  on  finite-difference  SGT  database 

In  the  SGT-based  approach,  the  displacement  can  be  expressed  as  (Zhao  et  al.,  2006): 

u„(r;r,)  - 

-|KG/.(r,;r)  +  ^G,(r,;r)]M„  (1) 

where  V  is  the  location  vector,  r,  is  the  source  location,  M,  i  the  moment  tensor  component  (ij=  1-3),  G(  r;r,) 

and  G(rj;r)  the  Green  tensors  for  a  unit  impulsive  force  at  the  source  and  T ,  respectively,  and  the 

strain  Green  tensor  composed  of  the  spatial  gradient  of  the  Green  tensors.  The  subscript  n  denotes  the  component  of 
the  displacement.  The  source-receiver  reciprocity  and  the  symmetry  of  the  moment  tensor  have  been  applied  in 

Equation  (1 ).  Notice  the  SGT  can  be  at  any  location  in  the  medium  and  thus  can  also  be  at  the  receiver  location  Tr  . 
Equation  (1)  is  the  basis  upon  which  the  FDSGT  moment  tensor  inversions  are  developed  (Zhao  et  al.,  2006;  Shen  et 
al.,  2009). 
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The  waveform  change  at  rR  due  to  a  local  model  perturbation  at  an  arbitrary  position  r  ,  can  be  expressed  as  (Zhao 
et  al.,  2006): 


<5w„(i*;r,)  =  HIJn(r;  rR  )dCiJkfj)elk  (rjr,)  (2) 

where  dCijfX r)  is  the  perturbation  of  the  elasticity  tensor  and  is  spatial  gradient  elements  of  the 

displacement  vector,  which  requires  one  forward  wave  simulation  from  the  source.  The  total  waveform  ehange  due 
to  the  model  perturbation  in  the  entire  volume  ean  be  found  by  volume  integration  of  Equation  (2)  in  the  model.  The 
above  waveform  perturbation  equation  is  the  essential  ingredient  in  the  SGT-based  tomographie  method.  The 
Frechet  kernels  relating  elastic  and  anelastie  properties  to  any  seismie  data,  such  as  frequency-dependent  traveltinic 
and  amplitude  anomalies,  ean  be  obtained  from  it  (Zhao  et  al.,  2005;  Zhang  et  al.,  2007). 

Figure  1  is  a  schematic  illustration  of  the  unified  tomography  and  source  moment  tensor  inversion  based  on  3D 
FDSGT  databases.  Station  SGT  databases  are  constructed  from  3D  reference  models  by  finite-difference  simulation 
of  the  responses  to  orthogonal  unit  impulsive  point  forces  acting  at  the  stations.  The  SGTs  in  a  small  volume 
surrounding  the  source  reference  location  extracted  from  the  databases  make  it  possible  to  invert  for  souree  moment 
tensors  and  location  in  a  global  optimization  seheme.  Forward  wave  propagation  simulation  from  the  souree  in  the 
3D  reference  model  provides  synthetic  waveforms  at  the  station  and  within  the  entire  model.  Travel  time,  amplitude, 
or  waveform  anomalies  are  measured  from  the  observed  and  synthetic  waveforms  at  stations.  The  forward  wave 
field  (from  the  souree)  and  the  station  SGTs  are  used  to  calculate  the  finite-frequency  sensitivities  to  perturbations  in 
Vp,  Vs  (or  bulk  and  shear  moduli),  density,  and  attenuation.  Together  with  the  sensitivities  to  souree  parameters,  the 
measurements  and  structural  sensitivity  kernels  are  used  to  invert  for  the  earth  structure.  The  inversion  results  are 
added  to  the  3D  reference  model.  This  process  ean  be  repeated  to  progressively  improve  the  resolution. 

Source  characterization  and  construction  of  hierarchical  FDSGT  databases 

Most  earthquake  souree  studies  to  date  utilize  1 D  earth  models  in  the  calculation  of  synthetic  waveforms  or  Green's 
functions.  To  accommodate  3D  structural  heterogeneities,  segments  of  body  and  surface  waves  are  time-shifted 
before  waveform  fitting  (e.g.,  Zhao  and  Helmberger,  1994;  Zhu  and  Helmberger,  1996).  In  places  with  highly 
heterogeneous  structure,  this  “eut-and-paste”  of  waveform  segments  may  not  fully  account  for  the  complexity  of 
wave  phenomena  caused  by  3D  structures  (e.g.,  wave  foeusing/ defoeusing).  One  of  the  major  sources  of  errors  in 
moment  tensor  inversion  is  phase  skipping  between  the  observed  and  synthetic  waveforms.  When  traveltime  delays 
due  to  3D  velocity  heterogeneities  are  comparable  to  the  wave  period,  phase  skipping  becomes  a  challenging 
problem. 

Advance  in  computation  has  made  the  use  of  3D  Green's  functions  in  earthquake  source  inversion  an  increasingly 
attractive  option  in  routine  operation.  Liu  et  al.  (2004)  determined  the  source  mechanisms  of  small  to  moderate 
earthquakes  in  southern  California  using  synthetic  waveforms  calculated  with  the  spectral-element  method.  Because 
the  deriv  atives  of  the  souree  parameters  are  determined  numerically  by  differentiating  synthetics  with  respect  to  the 
souree  parameters,  up  to  10  forward  simulations  are  needed  for  each  earthquake  (six  moment  tensor  components, 
latitude,  longitude,  depth,  and  the  reference  location).  When  the  reference  location  is  far  from  the  true  location,  the 
deriv  atives  of  the  location  parameters  may  not  adequately  reflect  the  non-linear  variation  of  waveforms  as  a  function 
of  location.  Taking  advantage  of  the  source-receiver  reciprocity,  Zhao  et  al.  (2006)  introduced  the  use  of  strain 
Green's  tensor  in  source  inversion  based  on  3D  reference  models.  We  have  extended  the  SGT-database  approach  to 
regional  and  teleseismie  scales  (Slien  et  al.,  2009),  so  we  are  able  to  use  any  combination  of  local  and  teleseismie 
waveforms  in  souree  inversion.  We  fit  filtered  and  time-shifted  waveforms  directly.  For  a  given  souree  location,  the 
waveforms  are  linear  functions  of  the  moment  tensor  elements,  w  hich  ean  be  solved  as  a  linear  inverse  problem.  To 
account  for  perturbations  in  earthquake  souree  location,  we  carry  out  a  grid  search  in  the  vicinity  of  the  reference 
location. 

Figure  2  shows  the  selected  earthquakes  in  the  southeast  Tibetan  plateau  and  the  stations  used  in  moment  tensor 
inversion  (Shen  et  al.,  2009).  We  calculate  synthetic  waveforms  using  a  non-staggered  grid  finite-difference  code  in 
a  polar-spherical  coordinate  system  (Zhang  et  al.,  2009).  Wave  simulations  are  carried  out  on  three  geographical 
scales:  loeal/near  regional,  far  regional,  and  teleseismie.  An  updated  crust  and  upper  mantle  model  (CUB  2.0, 
Shapiro  and  Ritzwoller,  2002)  is  used  as  the  3D  reference  model.  The  FDSGT  moment  tensors  agree  well  with  the 
global  centroid  moment  tensor  (GCMT)  solutions  (Figure  3),  with  one  notable  exception  at  the  eastern  Himalayan 
syntaxis  (event  03.230  in  Figure  3),  The  GCxMT  solution  shows  a  strike-slip  event  at  33  km  depth,  while  our  FDSGT 
solution  is  a  shallow  (2±  1 6  km)  thrust  event. 
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Figure  4  shows  an  example  of  the  fit  between  the  observed  waveforms  and  the  synthetics  calculated  from  the 
FDSGT  databases.  An  inspection  of  the  waveforms  in  detail  (Shen  et  al.,  2009)  shows  that  the  regional  stations,  the 
3D  model,  and  shorter-period  body  waves  than  those  used  in  the  GCMT  solutions  tend  to  provide  tighter  constraints 
on  the  source  location,  particularly  the  depth. 

Full-wave  tomography  and  finite-frequency  kernels  based  on  3D  reference  models 

To  date  the  inversions  of  P-  and  5-wave  speeds  have  been  carried  out  separately  under  the  assumption  that  P 
traveltimcs  arc  affected  only  by  the  P- wave  speed  of  the  clastic  media  and  5  traveltimes  by  the  5- wave  speed.  The 
practice  stems  from  ray  theory,  in  which  P  and  S  waves  arc  completely  decoupled  except  for  conversions  at 
discontinuities  following  the  Snell’s  law.  Based  upon  numerical  and  analytical  solutions,  we  have  shown  that  for 
Finite-frequency  waves,  5-speed  perturbation  may  have  significant  effects  on  P  waveforms  (Zhang  and  Shen,  2008). 
At  plaees  near  the  souree  or  receiver,  the  cross-dependence  of  P  waveforms  on  5-speed  perturbation  is  substantial. 
For  example,  the  integrated  absolute  P-travcltimc  sensitivity  to  5-spccd  perturbation  is  31%  of  that  to  P-spccd 
perturbation  when  the  source-receiver  distance  is  seven  times  the  wavelength  (L/X=7,  where  Z  is  the  wavelength  and 
L  the  distance  to  the  receiver).  At  a  source-receiver  distance  of  100  wavelengths,  the  ratio  is  about  5%.  This  Finite- 
frequency  effect  has  important  implications  for  seismic  tomography,  particularly  at  regional  and  local  scales,  where 
target  structures  are  often  located  w  ithin  several  tens  of  wavelengths  from  receivers  or  from  both  sources  and 
receivers.  Even  for  regions  far  from  sources  and  receivers,  neglecting  this  cross-dependence  means  that  the 
estimates  of  the  P  speed  near  the  souree  and  receiver  arc  likely  biased  and  such  errors  may  propagate  into  the  rest  of 
the  model  in  tomographic  inversion. 

Figure  5  shows  the  phase  sensitivities  of  a  regional  Rayleigh  wave  to  5-  and  P-wave  speeds.  Similar  to  the  results 
for  southern  California  in  Liu  and  Tromp  (2006),  there  arc  large  3D  variations  in  the  finite-frequency  kernels  due  to 
wave-speed  variations.  Although  commonly  neglected  in  surface  wave  tomography,  P- wave  speed  (and  density) 
affects  Rayleigh  waves  (also  sec  Zhou  et  al.  (2004)  and  Liu  and  Tromp  (2006)).  The  maximum  sensitivity  to  P- wave 
speed  is  on  the  surface  and  at  a  level  comparable  to  the  Vs  kernel.  It  decreases  to  -10%  of  the  maximum  at  depth  of 
~X/4.  Since  the  peak  of  the  Vs  sensitivity  is  at  the  depth  of  — X/3  and  heterogeneities  in  the  mantle  arc  generally 
much  weaker  than  in  the  crust,  it  may  be  necessary  to  include  the  sensitivities  to  P- speed  and  density  for  long-period 
(-50+  s)  Rayleigh  waves  in  a  high-resolution  inversion.  When  integrated  with  Pn ,  Sn ,  and  Love  waves,  Rayleigh 
waves  may  provide  additional  constraints  on  the  P- wave  structure  in  the  crust.  At  the  very  least,  incorporating  the 
effects  of  P- speed  on  Rayleigh  waves  may  minimize  a  likely  cause  contributing  to  the  apparent  Lovc-Raylcigh 
discrepancy  and  radial  anisotropy. 

Figure  6a  and  6b  show  the  P  and  5  wave  velocities,  respectively,  of  the  shallow  crust  of  the  southeast  Tibetan 
plateau.  The  model  is  derived  from  a  joint  inversion  of -10,000  measurements  of  phase  anomalies  of  P/Pn ,  S/Sn  and 
Raylcigh/Lovc  waves  in  three  frequency  bands  ranging  from  0.008-0.1  Hz  (Shen  et  al.,  2008b).  After  two  iterations, 
the  combined  velocity  anomalies  relative  to  the  initial  reference  model  (CUB  2.0,  Shapiro  and  Ritzwoller,  2002) 
show  that  the  observ  ations  require  P  velocity  increases  just  north  of  the  Indus  Yalu  suture  (IYS),  near  the  Bangong- 
Nujiang  suture  (BNS),  and  west  of  the  Xianshuhe  Fault  system  (XSHF).  These  high-velocity  anomalies  correlate 
well  with  areas  of  granitic  and/or  more  metamorphosed  Proterozoic  and  Paleozoic  rocks  among  the  regional 
Mesozoic  background  (Geological  map  of  the  Tibetan  plateau,  China  Academy  of  Geological  Sciences,  Chengdu). 
The  P-wave  resolution  at  such  a  shallow  depth  is  attributed  in  part  to  the  sensitivity  of  Rayleigh  waves  to  P-wave 
speed  perturbations  (Figure  5).  The  updated  5-vclocity  model  (Figure  6b)  shows  a  sharp  velocity  contrast  that 
follows  the  boundary  of  the  Tibetan  plateau  as  outlined  by  the  1000  m  elevation  contour. 

CONCLUSIONS  AND  RECOMMENDATIONS 

We  have  developed  a  unified  tomography  and  earthquake  moment-tensor  inversion  based  on  3D  reference  models. 
The  new  approach  accounts  for  complex  wave  propagation,  enables  fuller  utilization  of  an  arrival  on  all  three 
components  of  seismic  records,  and  linearizes  the  inverse  problem  by  iteratively  updating  the  3D  reference  model. 
We  calculate  synthetic  waveforms  using  a  non-staggered  grid  finite-difference  code  in  a  polar-spherical  coordinate 
system.  The  new  moment  tensor  solutions  accommodate  the  3D  velocity  heterogeneities  and  are  an  integral  part  of 
the  unified  tomography  and  source  inversion.  Body  and  surface  wave  traveltimcs  are  inverted  jointly  for  P  and  5 
wave  speed  perturbations.  The  results  show  a  much-improved  resolution  compared  to  the  reference  model  and 
regional  models  of  similar  scale. 
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Figure  1.  Schematic  illustration  of  the  unified  tomography  and  moment  tensor  inversions  based  on  3D 
FDSGT  databases. 
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Figure  2.  Selected  earthquakes  (red  circles)  in  the  southeast  Tibetan  plateau  during  2001-2004  and  the 
stations  (black  triangles)  used  in  moment  tensor  inversions.  White  rectangles  outline  the  three 
geographic  scales  of  the  finite-difference  nave  simulation:  (1)  local/near  regional,  (2)  far-regional, 
and  (3)  telescismie. 
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Figure  3.  Comparison  of  the  FDSGT  (white-and-blue  symbols)  and  GCMT  (whitc-and-red  symbols)  solutions 
for  13  earthquakes  during  2001-2004  in  the  southeast  Tibetan  plateau.  The  numbers  above  the 
FDSGT  solution  denote  the  y  ear  and  Julian  day  of  the  earthquake.  Also  plotted  arc  the  source 
mechanisms  of  10  additional  earthquakes  without  GCMT  solutions  (black-and-white  sy  mbols). 
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Figure  4.  (a-d)  Observed  waveforms  (thick  gray  lines)  are  compared  with  the  synthetics  for  the  best  source 
moment  tensor  solution  (red  lines).  The  frequency  bands  are,  from  left  to  right,  0.00833-0.0167, 
0.0167-0.0333,  0.0333-0.0667,  and  0.0667-0.167  Hz.  Dark  thick  lines  arc  time-shifted  observed 
waveforms  used  in  the  final  inversion.  For  the  first  arrivals  (in  c  and  d),  the  time  is  redueed  to  align 
the  arrivals  within  the  time  window.  The  amplitude  is  normalized  by  individual  trace. 

(e)  Normalized  w  aveform  misfit  as  a  function  of  depth.  Crosses  mark  the  grid-search  depths.  The 
shaded  region  indicates  the  la  confidence  limit.  (0  Contours  of  various  confidence  limits  at  the 
depth  of  the  global  best  solution  (8  km  for  this  event).  The  cross  marks  the  FDSGT  location;  the 
square  is  the  EHB  location  (Engdahl,  see  acknowledgements);  and  the  diamond  is  the  global 
centroid  moment  tensor  location.  From  the  innermost  contour  surrounding  the  FDSGT  solution,  the 
contours  are  the  68%  (dark  line),  80%,  90%,  95%  (dark  line),  and  99%  confidence  limits. 

(g)  Triangles  and  cross  mark  the  stations  used  in  the  source  inversion  and  the  event,  respectively. 
Gray  contours  follow  0,  2000  and  4000  m  elevation. 
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Figure  5.  Phase  sensitivity  kernels  of  a  regional  Rayleigh  wave  to  S  wav  e-speed  (top)  and  P  wave-speed 

(bottom)  perturbations  relative  to  a  3D  reference  model.  The  period  of  the  wave  is  10-20  s.  The  two 
vertical  eross-seetions  eut  the  model  along  the  great  circle  path  between  the  earthquake  source 
(in  the  southeast  corner  of  the  model)  and  receiver  (in  the  northwest  corner  of  the  model)  and 
perpendicular  to  the  great  circle.  Also  shown  is  a  quadrant  of  a  horizontal  section  on  the  surface. 
The  scale  is  ±2x10  13  s/m3. 
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Figure  6a.  P-wave  velocity  variations  at  6  km  depth  in  the  southeast  Tibetan  plateau.  The  initial  3D  reference 
model  (top  right)  is  converted  from  CUB  2.0  by  a  VplVs  ratio  of  1.74  in  the  crust  and  the  VpiVs  of 
AK135  in  the  mantle.  The  unit  of  the  color  bar  is  km/s.  The  left  panels  show  the  velocity 
perturbations  after  the  l*1  and  2nd  iteration  and  the  total  velocity  change  in  percent.  The 
middle-right  and  bottom-right  panels  show  the  updated  models  after  the  Ist  and  2nd  iterations.  The 
main  tectonic  features,  the  Indus-Yalu  suture  (IYS),  Bangong-Nujiang  suture  (BNS),  and  the 
Xianshuihei  Fault  system  (XS11F),  are  marked  in  the  upper  right  panel.  The  1000  m  elevation 
contour  outlines  the  southern  margin  of  the  Tibetan  plateau  and  the  boundary  of  the  Sichuan 
basin. 
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Figure  6b.  Same  as  Figure  6a,  except  showing  S  velocity. 
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